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Abstract. The continuum limit of a recently- proposed model for charge transport in resonant- 
tunneling semiconductor superlattices is analyzed. It is described by a nonlinear hyperbolic inte- 
grodifferential equation on a one-dimensional spatial support, supplemented by shock and entropy 
conditions. For appropriate parameter values, a time-periodic solution is found in numerical simula- 
tions of the model. An asymptotic theory shows that the time-periodic solution is due to recycling 
and motion of shock waves representing domain walls connecting regions of the superlattice where 
the electric field is almost uniform. 
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1. Introduction. Self-sustained oscillations of the current have been observed 
in recent experiments on n-doped weakly-coupled semiconductor superlattices (SLs) 
under dc- voltage bias along the growth direction fl6|| . The main charge transport 
mechanism is quantum resonant tunneling through the SL [pd| |l7| , fL2|| , Similar but 
damped oscillations have been reported earlier in undoped superlattices subject to 
laser photoexcitation J2(J. In both cases a phenomenological discrete drift model 
proposed by one of the authors |6) explains the oscillations in terms of the formation 
and dynamics of electric field domains, i.e. regions in space in which the strength 
of the electric field is nearly uniform. These phenomena are thus examples of time- 
dependent pattern formation in resonant-tunneling SL, a scarcely explored area of 
semiconductor physics. 

We shall study asymptotically and numerically the self-sustained oscillatory solu- 
tions of the continuum limit of the model that applies when the number of quantum 
wells in the SL is large or equivalently when the SL is long. The continuum limit of 
the model, to be derived below, is given by the equations (see [Q) 

(1.2) y f E(x,t)dx = <f>. 

L Jo 

The unknowns in these equations are the current density I(t) and the electric field 
E(x,t) on a one-dimensional SL < x < L. The velocity v(E), shown in the inset 
of Figure la, is a positive function for E > having two peaks at E = 1 [v(l) = 1] 
and E = Em > 1 [v{Em) = vm > 1] separated by a minimum at E — E rn G (1,Em) 



♦Received by the editors of SIAM J. Appl. Math, on 11 July, 1995. Manuscript number 028888-1. 

^Escuela Politecnica Superior, Universidad Carlos III de Madrid, Butarque 15, 28911 Leganes, 
Spain. The research of this author was supported by DGICYT grants PB92-0248 and PB94-0375, 
by North Atlantic Treaty Organization traveling grant CRG-900284, and by EEC Human Capital 
and Mobility Program grant ERBCHRXCT930413. 

* Department of Mathematics, Duke University, Durham, NC 27708, USA. The research of this 
author was supported by the Army Research Office grant DAAH04-93-G-0011, National Science 
Foundation grant DMS-91-03386, and by North Atlantic Treaty Organization traveling grant CRG- 
900284. 



1 



2 



L. L. BONILLA, M. KINDELAN, M. MOSCOSO AND S. VENAKIDES 



[v(E m ) = v m < 1]. The constant ^ is a control parameter proportional to the dc 
voltage bias. The self-sustained oscillations are stable time-periodic solutions that 
appear for 1 < <f> < E m ; see the oscillation of the current density in the inset of 



Figure lb. Equations ( |l.l|) and ( |l.2|) are to be solved with the boundary condition 



(1.3) __ = c>0 . 

Eq. ( tl.lj ) is hyperbolic, so its solutions will in general develop shock waves. We 
supplement the model with the shock condition |ij: 
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and the entropy condition 

(1.5) v{EJ)>V{E+,EJ)>v{E+). 

Here V(E + , E_) is the velocity of a shock wave located at x — X(t) and such that 
E(X-,t) = E_ and E(X+,t) = E+. 



Physically, Eq. (L2) expresses the fact that the volt age drop across the semicon- 



ductor length is kept constant at the value 0L. Eq. (1.1) can be derived directly from 
Ampere's Law 

OE 

(1.6) —+v(E)n = I, 

by substituting in it the Poisson equation 

(1.7) 

where the one in the right hand side represents the normalized density of the dopant. 
The electron current v(E) n is essentially given by the probability of tunneling from 
one quantum well of the SL to the next one, and is assumed to be proportional to 
the electron density in the quantum well. Our derivation of fll.lD homogenizes over 
the SL structure. The quantum effects due to the SL structure that dominates the 
phenomenon are modeled in the shape of the function v{E) (the electron velocity) and 
in particular in its having two peaks as shown in the inset of Fig. la above. Peaks 
occur at certain values of the field as a result of energy level alignment in adjacent 
quantum wells that leads to enhanced (or resonant) tunelling and hence to enhanced 
transmission. 

The shock condition ( |1 .4[ ) arises in the passage from the discrete model to the 
continuum model. It has the geometric interpretation of an equal area rule as one 
observes in Figure 2. In essence the shock condition provides a relation between the 
field E+ at the front of the shock wave, the field E- at the back of the shock and the 
shock speed V . The domain of this relation is restricted by the requirement that the 



entropy condition be satisfied. The entropy condition (1.5) is the usual one, saying 
that the particles entering the shock move faster than the shock itself, which in turn 
moves faster than the particles it encounters. The entropy condition also ensures that 
the electron density inside the shock wave (see below for the meaning of a region of 
the SL inside of the shock wave) is positive. 
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Our calculation is organized as follows: In Section || we discuss the oscillations 
from a qualitative viewpoint and in the further sections we do the analysis of one 
complete cycle of the oscillations. We start by analyzing the Shock Kinematics, in 
Section [|, and explain analytically the numerically observed monopole (shock plus 
tail) structure. In section |] we introduce our asymptotic scaling, derive the outer 
and inner solutions and put the pieces together working and matching in sequence 
the different regions which describe a whole oscillation cycle. In sections || and || we 
analyze in detail the process of creation and disappearance of a monopole, and in 
section ^ we include some final remarks. Finally Appendix [A| contains a derivation of 
the continuum limit of the discrete SL model. 

2. Qualitative Discussion of Current Oscillations. We study the asymp- 
totic limit in which the length of the SL sample is large. We have derived the time- 
periodic solution of the above model asymptotically and numerically. We describe it 



in Figure 1, in which the numerical solution of equations (1.1)- (1.3) is represented. 
This numerical solution has been obtained by the method of lines, using backward 
differences in space, a fourth order Runge-Kutta method in time, and a large number 
of grid points (5000) in order to have sufficient precision in the numerical approxi- 
mati o n. No tice that this method is, in fact, equivalent to solving the discr ete m ode l 



( A.7)-( A.1C ) as an initial value problem. We have also solved equations (fLl])-( |L~3| ) 



using a fully implicit finite difference scheme and including an artificial diffusion term 
proportional to v(E). The key element of this method is a modification of the method 
proposed in which distributes the grid points according to a monitor function sen- 
sitive to the local values of the slope of the field distribution, and tuned to track the 
travelling internal layers. The results obtained with both methods agree. 

One observes that at all times the field is essentially uniform over two or three 
spatial intervals called domains. The latter are separated by one or respectively two 
thin internal layers (IL) in which the field value experiences a positive jump. The 
IL(s) are monotonically increasing in x. To follow a complete cycle, we start with two 
domains separated by a single IL. Before the IL has reached the end of the sample 
second IL starts forming at an interior point to the left of the first one, and 
also moves to the right. The field now displays two ILs separating three domains. As 
time still increases, the IL on the right disappears either because it reaches the right 
end x = L of the sample (for 1 < <f> < cf>d), or because the jump between the field 
values to its left and to its right decreases to zero (for <fi > (j>d, where 1 < (j>d < E m ). 
This leaves a field profile that has only two domains. The cycle then repeats itself. 
Throughout the cycle, the uniform value of the field in each domain varies continuously 
and periodically with time. At each time, the field is an increasing function of the 
spatial variable. See Figure 1. 

In each IL, we find that a shock connects the left or right value of the field 
with an intermediate field value. The latter, in turn, connects to the value on the 
other side of the IL by an exponential tail that moves rigidly with the shock. The 
characteristics are parallel to the shock line on the side of the tail and end into the 
shock line transversely on the other side. We call the combination af a shock and a 
tail a "monopole wavefront" . In this language, an internal layer (IL) is a monopole 
wavefront that consists of a shock and a right tail or a shock and a left tail. As we 
explain in section ||, the structure of a shock with a left AND right tail does not 
appear. 

As the field displays the time-periodic spatial pattern formation described above, 
the current performs a simple oscillation as seen in the inset of Figure lb. Roughly 
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speaking, it increases with time as long as there are only two field domains and 
decreases when there are three field domains. We can describe how the field values 
on the domains depend on the current without reference to the position of the shocks 



in the following way. In the asymptotic scaling of interest [see Eqs. (4.1) to (4.4) 
below] , in which the length of the semiconductor sample and the speed of the shocks 
are of order 0(1), the quasistatic approximation / = v(E) is valid as long as the 
current value is between the local maximum and the local minimum of the v curve 
(see the inset of Fig. la) and not too close to these extreme values. The three roots 
of the equation above are denoted by E^'(I) < E^(I) < E^(I). As the current 
increases [region I, times (1) and (2) in Figures 1] taking values in the interval (v m , 1), 
there are only two domains separated by an IL. The field on the left domain takes 
the value E^(I), while on the right domain it takes the value E^(I). The shock 
condition and the entropy condition, plus the requirement that monopoles can have 
only a left OR (exclusive) right tail allow us to determine uniquely the shock and 
tail characteristics as functions of the current. After the current reaches its maximal 
value, a value that, as we will see, is a little higher than the peak of the v curve, a new 
shock is created (region II, time (3) in Figure 1). As the current value now decreases 
(Region III, time (4) in Figure 1), there are three field domains. The field value on 
the left domain is E^(I), the value on the middle domain is E^ 2 \l), and on the right 
domain it is E^ (I). Again, using the shock and entropy conditions, we can determine 
uniquely the characteristics of the two shocks and tails. This determination is shown 
graphically in the five pictures in Figure 3 that describes the whole cycle. We are 
forced to graph 1/v instead of the more natural v in this figure, because v appears in 
the denominator of the shock condition. Thus, for example, instead of picturing the 
current increasing we think of 1 // decreasing. The monopole on the right disappears 
(Region IV; not shown in the figures) as the current is close to its minimal value, a 
value that is equal or slightly higher than the minimum of the v curve, depending on 
the bias. 

The position of the IL in region II is directly determined from the voltage bias 
condition as a function of /. Since the shock position and speed are known functions of 
the current, we easily obtain an ode for the current as a function of time from which 
the time-dependence of the current and of the shock position may be derived. A 
similar but more involved calculation holds in region III, where the current decreases 
taking values in the interval (v m ,l). We have three relations: a formula for each 
of the two shock speeds as functions of the current and the voltage bias condition. 
The latter is a linear relation between the positions of the two shocks in which the 
coefficients are known functions of the current. We have three unknowns, the two 
IL positions and the current. We can determine the unknowns as functions of time 
given appropriate initial conditions. We provide the required initial conditions by 
(a) calculating the time and position of shock formation and (b) matching with the 
field in region II, in which the current is near its peak value. In regions II and IV 
the quasistationary approximation is not always valid. We use fast scale variables to 
describe the birth and death of a shock shock in detail and match them to the regions 
I and III. 

We can represent the evolution of the shock wave without consideration to its 
spatial trajectory on the SL by means of the following diagram. We define the function 



E + = ^(E-tu) (where u = 1/V) by solving the equal area rule (L4) for E + . This 
function yields the field in front of the shock wave, given appropriate values of the 
shock speed and of the field at the back of the shock. The domain of this function 
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consists of all the points (E_, 1/V) for which there is an E + such that the equal area 
rule and the entropy inequalities are true. Thus a point in the domain characterizes a 
shock completely. In Fig. 4 we have depicted the domain of T(E- , u) and the complete 
life-cycle of a shock wave in it. As one can see in the figure and easily derive, the 
domain of !F(E-,u) is bounded below by the graph of the function l/v(E). It is 
bounded above by a line on which E+ reaches the maximal value obtained by the 
equal area rule that does not violate the entropy inequalities: = V(E+, E-). 

This line intersects the graph of l/v(E) at point (Ep, 1/-Tf) on the first branch. It 
again intersects the graph at the common point of the second and third branches 
(E m , l/v m ). Observe that on the part of the first branch that belongs to the domain, 
we have v(E_) = V(E + , E_), while on the second branch we have E + = E_. Also 
observe that E_ = E^\j F ), and E + = E^(I F ). 

According to both our numerical experiments and asymptotic calculation, the 
shock wave moves so that the point (E-, 1/V) always remains on the boundary of the 
domain of F(E-,u). The shock wave first appears at a point on the second branch 
of l/v(E) where E+ = E- (see Figure 4). Then the point (E— , 1/V) moves from the 
second to the first branch of l/v(E), where it satisfies V(E+, E-) = v(E-) until it 
arrives at (E F , I/If)- It then continues moving on the upper boundary of the domain 
of !F(E-,u) where V(E +> E-) = v(E + ), until it reaches (E m , l/v m ) where the shock 
wave dies. (We have assumed a large enough bias, <f) > <f)d- For smaller biases, the 
shock reaches x = L before the current has attained its minimum value. Then the 
shock dies at some point of the upper boundary of the domain of !F(E-,u)). 

To summarize, in a full cycle of the oscillation the following regions can be iden- 
tified: 

• Region I: One monopole (two domains) exists and the current is increasing 
(see Figure 3.e). 

• Region II: A second monopole is born near x = and the current is nearing 
a maximum which is lightly higher than the peak of the v curve (see Figure 
3i). 

• Region III: Two monopoles (three domains) propagate and the current de- 
creases (see Figure 3.b). 

• Region IV: The monopole on the right disappears and the current is nearing 
a minimum which is slightly higher than the minimum of the v curve (see 
Figure 3.c). 

The different regions here described are analyzed asymptotically in section [|, 
where we derive the outer and inner solutions and put the pieces together working 
and matching in sequence on regions I, II, III and IV. The technical work on regions II 
and IV is done in sections || and ^ respectively. 

3. Kinematics of Shock Waves . In this section we shall analyze the pro- 
cesses of shock formation (near the boundary) and the existence and properties of 
monopole wavefronts which are a shock wave and a tail to its left or right moving at 
the shock velocity. These matters are fundamental ingredients of the asymptotics of 
later sections. 

3.1. Formation of Shock Waves . The formation of shock waves on the infinite 
real line from arbitrary initial data and constant current have been considered before 
for the model of the Gunn effect in semiconductors Jl8|, |2[ ^] . The main result is that 
for an initial positive field profile a shock wave develops in finite time if dE/dx > 
is large enough. For the time-periodic solution we are interested in, the shock waves 
develop near the boundary at x = for appropriate values of the current density. We 
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shall now analyze Region II and derive exact formulas for the time and position at 
which a new shock wave appears. 



Let us solve the hyperbolic equation (1.1) with the boundary condition (|l.3| ) by 
the method of characteristics. 



(3.1) 
(3.2) 

with the boundary conditions: 

(3.3) 
(3.4) 



dE(t; t) 
Jt 



= I-v(E(t;r)), 



dx(t; t) 
Jt 



v(E(t;r)), 



E(t;t) =E (t), 
x(t; t) = 0. 



We have parametrized the characteristic curves by the times at which they issue 
from the boundary at x = 0. The boundary field E (t) has to be chosen so that the 
boundary condition dE(0,t)/dx = c holds. Clearly, 



(3.5) 



dE .„ , BE , . 8t , 
c =^zM = —(r;r)—\ x ^ 



dx 



dr 



dx 



We derive some relations that will be of use. Taking a r-derivative of the solution of 
Eqs. @ and Q 



(3.6) 



v (E(s; r)) ds, 



and then setting x = (equivalently, t = r) in the result, we obtain 



(3.7) 



^■(t;t) = -v(E(t;t))=-v(E q (t)) 



Eqs. (3.5) and (3.7) yield 



(3- 



dE 



(r;r) = -cv(E (t)). 



An equation for Eq is found by taking a r-derivative of the boundary condition ( |3.3| ) 
and substituting (3.8) in the result: 



dE 
dr 



dE BE 

-fo^ T ) + a7( r ' T ) = J ( T ) - "(^W) - cw (£o(t)), 



that is 
(3.9) 



dE Q ( T ) 
dr 



+ (c+1)v(E (t))=I(t). 



This equation may be obtained directly by substituting c = dE /dx into the hyperbolic 
equation (PI). Finally, we can combine ( |3.l| ) and ( |3.2| ) into the equation 



(3.10) 



x{t; t) + E{t; r) = / I(s) ds + E (t) 
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We can now derive a formula for the time at which a shock wave first appears. This 
happens when the function r(x;t), obtained by solving x(t;r) = x for fixed t, first 
becomes multivalued. At the time of shock formation, the function x(t;r), has an 
inflection point with zero slope in the r variable and therefore satisfies dx(t; t)/8t = 0, 
and d 2 x(t;T)/dr 2 = 0. These are two conditions to be used in determining the time 
and location of the new shock wave. 



(3.11) 



dx(t; t) 
dr 



dE dE(t;r) dE(t;r) 

w Kr) = -(1 + c) v(E (r)) - 



where£3.9) has been used. We can find dE(t; t)/8t = E T (t; r) by differentiating (3.1) 
and (3^) with respect to r, using (:?■£), and solving the resulting linear problem: 



(3.12) 
(3.13) 



dE T (t;r) 
dt 



= -v' {E(t; t)) E r (t; T ), 
E t (t;t)=-cv{E (t)). 



The result may be inserted into (3.11) to obtain 
dx(t; t 



(3.14) 



v(E (t)) c exp 



v'(E(s;r))ds 



From Eq. (3.14) we immediately obtain 

d 2 x(t;r) v'(E (r)) _ , ,dx(t;r) 



dr 2 



v(E (r)) 



(3.15) x[v'(E (t))+cv(E (t)) / v"{E{r,T)) exp 



Ot 
t 



+ cv(E (t)) exp 



v'(E{s;T))ds 



v'(E(s;T))ds 



dr 



When the expressions given by Eqs. ( 3.14 ) and ( 3.15 ) are set equal to zero, their 
solutions t = t s and r = r s yield the time and position at which a shock wave first 
appear, 



(3.16) 
(3.17) 
(3.18) 



1 + c 



exp 



v'(Eq{t 3 )) 
~cv(E q {t s )) 



v"(E{r,T s )) exp 



v' (E(t; Ts )) dt 
v'(E(s;T s ))ds 



dr. 



v(E(t;T s )) dt. 



The times t s and t s in Eqs. (3.16) and (3.18) are functions of c and L determined 
by the condition that t s is the smallest time for which (3.16) and (3.17) are satisfied. 
Notice that (3.16) has no solution t s > t s for — 1 < c < 0. This explains why the 
time-periodic electric field profiles, which appear in the numerical simulations with 
such boundary condi ti ons d o not exhibit shocks In order to find t s , x s and r s 

from equations ( 3.16 )-( 3.18 ), we need to know E(t; t) along the characteristics, which 
is obtaine d by solving (B.l) when the current is known. This analysis is carried out 
in Section 4.3 using an asymptotic approximation to the current density I{t). 
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3.2. Monopole Wavefronts. We now construct the ILs refered to in the In- 
troduction, which are basic elements of our asymptotics. Assuming constant current, 
we construct shock waves with a left or right tail moving rigidly at the shock velocity 
V(E+, E-). We call a monopole with a left (resp. right) tail the internal layer com- 
posed of the shock wave plus the transition layer (the tail) to its left (resp. right). 
The monopoles will be used as building blocks in the asymptotic description of the 
time-periodic solution. In the moving coordinate x = x — X(t), where X(t) is the 
position of the shock at time t [so dX/dt — V(E+, E—)], we have 

Fie 

(3.19) [v(E)-V(E + ,E-)] — =I-v(E). 

ox 

We shall distinguish two cases: 

Case 1: I < V. We construct a left tail by solving fl3.19| ) in the region x < 0. As x 
varies from to x — — oo, E varies from E(0—) = E_ to E{— oo) = El- The quantity 
v(E) varies from v(E-) to v(El) = I [El is necessarily a fixed point of ( 3.1 9| ) ] . Let 



us prove that u(£L) = V{E + ,EJ). If V < v(EJ), we have v(E L ) = I < V < v(E-), 
and the function v(E) — V has a zero between El and E_. Thus ( |3.19| ) does not have 
a solution that connects El and £L, and we must have v(E_) < V(E+, E-). But 
then the entropy condition V(E + ,EJ) < v{EJ) implies V (E + ,E-) = v(E-). The 
behavior of the tail near x = is described by 



dE I -V 



d x v>{E_){E-E-) 
which when solved yields 



E(Q)=E_, 



With this choice of the minus sign in front of the square root, E(x) is an increas 
ing function, as indicated by ( 3.1 9| ) . Since V > I, v'(E) > and therefore E. 



is 



necessarily on the first branch of the v curve. 



The only solution of (3.19) for x > with E(0+) = E + that we have observed in 
the numerical simulations is the constant solution E(x) = E + , with v(E+) — I. We 
have the following partial explanation of this fact. Let E (x) — > En as x ~~ * 00 ■ Clearly 
v(Er) = I. In a neighborhood of Er we can multiply (§T|) by v'(E) - v'(E R ) and 
obtain 

d , n v'(E R ) , n 

As x ~ * 00 7 the point v = I is unstable if v'(Er) > and stable if v'{Er) < 0. The 
unstable case corresponds to fields on the third branch of the v(E) curve, and the only 
solution having E — > Er [hence v(E) — ► v(Er)] as x — * oo is E = Er = E + (v = I). 
In the stable case, we cannot exclude right tails based on the above arguments. We 
have excluded them because they do not appear in our numerical simulations. 

Case 2: I > V. This case is completely analogous to case 1. We obtain monopole 
waves with a right tail. To summarize: 

• When I < V, u(£L) = V(E + ,E_) = U{E+), I = v{E + ), the monopole 

moves with velocity dX/dt = U(E + ), and it has a left tail, i.e., El < E- < 

E+ = Er. 
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When I > V, v{E + ) = V(E + ,E_ 
moves with velocity dX/dt — W(E_ 
E+ < Er. 



= W(E_), I = v(E_), the monopole 
and it has a right tail, i.e., El — < 



The functions U(E + ) and W(E_) are defined by the above relations. 

In Figure la a monopole with a right tail can be observed in the field profile 
corresponding to time (2). This is better displayed when plotting dE/dx as shown 
in the inset of Figure lb. Time (1) corresponds to a value of the current close to If 
(If = 0.58) where neither left nor right tail exists. 

4. Asymptotic Analysis of the Time-Periodic Solution. Let us describe 
one period of the oscillations for bias 1 < 4> < E^(l) (recall that E^(I) < E^ 2 \l) < 
E^ 3 >(I) are the three roots of v(E) — I for v m < I < 1). It is convenient to redefine 
the time and space scales in such a way that the SL length becomes 1: 



(4.1) 



1 



V= L> 



Then Eqs. (1.1) - (TLJ) become 



(4.2) 
(4.3) 
(4.4) 



dE 
97 



v(E) 



dE 
dy 



t 

S= L 



v(E) 



E(y,s)dy = <j>. 



8E(0, s) 
dy 



cL. 



We shall describe the time-periodic solution of these equations in the limit e — > 
(L — * oo) by leading-order matched asymptotic expansions. Observe that L is 
the nondimensional length of the superlattice [see Appendix |A|, equations ( ]A.6| ) and 
(A. 11)], and therefore we analyze the limit in which the SL length is large compared 
to e E M /(N D e). 

4.1. Outer solutions. Clearly the leading order of the outer expansion of the 

yields: 



solutions to Eqs. (4.2) 



(4.5) 



I -v(E) = 0. 



Then the outer electric field is a piecewise constant function whose profile is a succes- 
sion of zeroes of / — v(E), E^ k \l) (k — 1,2,3), separated by discontinuities. These 
discontinuities are the shock waves corresponding to monopoles with left or right tails 
moving with speeds given by (1.4), as discussed in Section [| Let us suppose that 
y = Y(s) represents the location of a shock wave separating two different solutions of 
(4.5). According to Section 0, 



= U{E+) 
W(E_), 



V(E + ,E_) 
V(E + ,E_) 



= v(E_) = U(E + ) if / < V(E + ,E_), 
v(E+) eeW(E-) if I>V(E + ,E_), 



Notice that the limiting values E + in (4J3) and £L in (4/7) are solutions of (45), and 
therefore they are functions of I(s). The cur rent density should then be determined 
from these equations and the bias condition (4.3), (see below). 
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4.2. Inner solutions. Near y = or near the shock waves there are regions of 
fast variations of the elec tric field. In them I = I(s), E ~ F(x, t), where x — y/e and 
t = s/e according to (4.1). 

The field in the boundary layer near y = obeys the equations: 



(4.8) 
(4.9) 



dF , x dF 

dEn 



F(0,t) =E (s), with e-^ + fl 

as 



I(s)-v(F), 
c)v(E ) = I(s), 



because of (|3.9| ). Except in very short time intervals when a new shock wave is 
being formed, F = F(x) is a quasi-stationary monotonically increasing profile joining 
F(0) = E^[I(s)/(l +c)] and F(oo) = E^[I(s)}. 



In the tail regions of a monopole, the electric field is a solution of ( 3.19 ), that is, 
F = F(x; s) with x = (y — Y(s))/e and 



(4.10) 



[v(F) - V(E, 



dx 



I(s)-v(F). 



When I(s) < V(E+,E-), [resp. I(s) > V(E + ,E_)), E_ (resp. E + ) is a function of 
E+ (resp. E-) given by V(E + ,E-) = v(E-) = U(E+) [resp. V(E+,E_) = v(E+) 



W(E-)], and ( 4.10| ) should be solved for x < (resp. x > 0) with the boundary 
condition F(0— ) = E- [resp. F(Q+) = E + }. Obviously F matches the value of the 
outer solution, E^ [I(s)], as we leave the tail region, x ~ * — 00 (k = 1,2) or x — > +oo 
(fc = 2,3). 

4.3. Putting the pieces together. We shall start at a time where there is only 
one shock wave on the SL, the current is If (see Fig. 4) and the electric field is a shock 
wave joining spatially uniform regions with E_ = EW(I F ), E+ = eW(I f ). At time 
s = 0, the shock wave is located at y = Y F = [E^(I F ) - (f)]/[E^\l F ) - E^(I F )] 
according to the bias condition (4.3). For s > 0, the shock (E-, l/V) is at the upper 
boundary of the domain of E + = T(E^,u) (see Introduction), so that v(E + ) = 
V(E+, E-) = W(E-), I < 1 and a monopole with a rigidly moving right tail is 
formed. Let y — Y(s) be the shock position. Then the outer field profile is 



(4.11) 



E(y,s)=EW (/( S )), 
E(y,s)=E^ (1(a)), 



if 
if 



0<y<Y(s), 
Y(s) <y<l. 



The bias condition ( f4~3| ) determines Y(s) as a function of I(s): 

(i)-4> 



(4.12) 



Y = 



E( 3 )(I)- EW(I)' 



Inserting (4.12) into (4.7) we find the following autonomous equation for the current 
density: 



(4.13) 



dl 



(E<® - EM) 3 v^W! 



> 0. 



ds E&) -4>+ ((f> -£■(!)) v'Jv's 
where w' = v'(E^), and Wj = W(E^), j = 1, 2, 3. To obtain (|4.13[) we have used 



dE^(I) 



ds 



1 dl 

— , 71=1,3, 



v'(E( n ) (I)) ds' 



Traveling fronts in semiconductor superlattices 



11 



which follows from time-differentiation of fl4-fi|) with E = E^ (I). Eq. ( 4.13 ) explicitly 
displays the quasi-st eady growth of the current during this stage of the oscillation. 
The solution I(s) of ( 1.13 ) reaches 1 at a time s = si which is numerically calculated. 
As w e approach this time, i? (1) (/) ~ 1 - [2 (1 - 7)/|u"(l)|] 1/2 , and the solution of 
(1.13) becomes 



(4.14) 
(4.15) 



1 - 



K(i)| 



[E^{1)-1] 2 W(1] 



£( 3 )(l)-0 

The corresponding outer electric field profile is 



(4.16) 



with 



(4.17) 



E(y, Sl ) = 1, 
E(y, Sl )=E^(l), 



if 
if 



0<y<F( Sl ), 
Y( Sl )<y<l, 



E^(l)-cf> 
EW(1) - 1' 



Obviously after this instant our approximations break down. What happens then? 

We shall see below that at s — sm, sm — si — 0(yfe), I reaches a maximum and 
then decreases, while the field to the left of the monopole, El — E-, increases linearly 
with time, and the field to the right of the monopole, E R = E^(I), decreases. When 
I surpasses 1, El can no longer be approximated by E^'(I). In a short time interval 
about s = sm, El and I are close to 1, and the difference I — v(El) eventually 
acquires a positive value of the same order as the time derivative edEL/ds = 0(e). 
Then I - 1 = 0(e), E L — 1 = 0(e 1 / 2 ), which happens in a time scale s—s M = 0(e 1 / 2 ). 
We thus make the ansatz 



(4.18) 

(4.19) 
(4.20) 



sm 



Er 



1 



,1/2 



El(S), 



~ e l/2 ' 

J~ l + e/(s), 
E R -E^(1) ~ e E R (s). 



In the time scale (4.18), the shock speed is 0(e 1 ^ 2 ) 

dY 



(4.21) 

so that 
(4.22) 



ds 



= e x ' 2 W(l)+0(e). 



Y(s) = Y(s M ) + e 1 ' 2 W(l) s + 0(e). 



[Notice that Y(sm) — Y(si) = 0(e)]. Inserting the field profile (4.20) and the time 
scale Q4.18D into Eq. (O), we get 



(4.23) 
(4.24) 



dE R 
ds 



I 



v' (EW(1)) E R . 
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We therefore have 
(4.25) 



Er 



to leading order. Inserting the electric field profile and (4.22) into the bias condition 
(4.3), we find E L by equating terms of order e 1 / 2 : 



(4.26) 



Er 



(E^(l)-l)W(l) 
Y( Sl ) 



a s, 



where a is given by (4.15). To find /, we substitute this result into (4.23), thereby 
getting 



(4.27) 

In outer units we therefore have 

(4.28) / ~ 1 + ea 

(4.29) 



f K(l)l 2-2 

1 = a a s . 



1 - 



\v"(l)\a(s-s M f 



E L 



2e 

1 + a (s - sm)- 



We have chosen E-(sm) = 1, as for this value of the electric field the current 
density reaches its maximum. The relation between sm and previous times for which 
the approximations ( 4.1l| )-( 4.13| ) hold should be found from a matching condition. 
I(s) takes on the value 1 — fte, K = O(l), at the times s K % and s K 2 1 given by 



(4.30) 



sm + (-iy 



l2(a + n)e 
a 2 \v"{l)\ 



with n = 1, 2. s\ for which I(s±) — 1 corresponds to n = in ( 4.30 ). At these times, 
E R - E(V(l) and 



(4.31) 



E L (s Kn ) ~ 1 + (-1)' 



'2(a + K)e 
W 



Notice that (4.28) and (4.14) match for any positive k = 0(1), as their difference is: 
(4.32) V2 |v"(l)| a 2 (a + k) e (s - Sl ) + 0(e) , 

which is 0(e 1 ^ 2 ) — o(l) when the proper outer scales of current [I = 0(1)} and time 
[s = O(L)] are used. We will choose appropriately the value of k (k = 1.78), so as to 
optimize the leading-order asymptotics we are describing here. 

During a time interval about s = Sm-i the time derivative of El is not negligible, 
and we have an unsteady stage during which a new shock wave is crea ted if the excess 
charge c is positive. Using the current density given by equation ( 4.2S ), the new shock 
wave appe ars at x = x s , t = t s (s s , y s in the outer variables) as obtained by solving 
Eqs. ( 3.16 ) - ( 3. IS ). Details about the birth of the new wave are given in Section [|. 
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After the new wave has appeared, there is a time interval where both old and new 
shocks (located at y = Y and y = Y n , respectively) coexist. The outer field profile 
consists of two shock waves connecting regions where the electric field is uniform: 



E = E Ll 
E = Em, 
E = E R , 



0<y<Y n , 
Y n <y<Y, 
Y < y < L . 



(4.33) 

The bias condition (JO]) is now 

(4.34) Y n E L + {Y -Y n )E M + {l-Y)E R 



Eq. (4.5) implies that 

(4.35) E L =EW(I(s)), E m =eW{I{s)), E r = E^{I(s)). 
The equations of motion for the monopoles are 

(4.36) 
(4.37) 



dY 

= U(E M ) 

as 



dY 

— = W{E M ). 
as 



From Eqs. (4.34) - (4.37), we determine the unknowns I(s), Y n (s) and Y(s) as 
follows. The bias condition ( ff.34 ) may be rewritten as 



(4.38) 
where 

(4.39) 



Y n = p-aY, 



or using (4.35), 
(4.40) 



E R — Em 
Em — El 



E (3) _ E (2) 



Er 



Em — El 



(3 = 



E^ - 
EW - E(V 



The three remaining unknowns Y, Y n and I can be determined from Eqs. (4.35) - 
(4.39). By using these equations, we can first find Y and Y n as functions of J, and 
then obtain an equation for Z(s). Eqs. ( 4.36 ) and ( 4.37 ) yield 



(4.41) 



dY n _ U 2 dY 
dl ~ W 2 dl ' 



A linear equation for Y(I) may be obtained from (4.38), (4.40) and (4.41): 

a' W 2 „ f3' W 2 



(4.42) 



dY 



Y 



U 2 + aW 2 ' U 2 + aW 2 ' 
where a 1 and /?' are the derivatives of a and [3 with respect to /. 
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A careful computation shows that the solution obeying Y(l) = [E^ (1)— <p\/[E^ (1)— | 
1] (leading-order location of the old shock wave at the time the new shock forms), 
and thus matches the previous stage, is 



Y = lim 

7^0+ 




(4.43) 



W 2 (3' 
U 2 + aW 2 

eW(i) -4> 

£7( 3 )(1) - 1 



exp 



exp 



a' Wo 



Un 



+ aW 2 
a'W 2 



■ dr 



dJ 



1-7 



U 2 + aW 2 



■ dr 



Eq. ( 4.43 ) holds as long as 

(4.44) < Y < 1 



and 



< I < 1. 



Taking a time derivative of the bias condition and using (4.36) and (4.37), we find 

dJ _ U 2 + aW 2 
ds 



(4.45) 



P'-a'Y ' 

to be solved with the initial condition I{s s ) = I s , (the value of the current at the 
shock formation time), which comes from matching with the previous stage. 

The solution of this equation shows that the current decreases until one of the 
two conditions (4.44) breaks down. Either (i) the shock wave at Y(s) reaches y = 1 
at some time Sd, or (ii) the current I(s) takes on the value v m corresponding to the 
minimum electron velocity at some time s m with Y(s m ) < 1. One possibility or the 
other is realized acco rding to the value of the bias: Let 4>d be the bias for which Y = 1 
when / = v m in Eq. ( 4.43| ). For 1 < <j) < <frd we have possibility (i), whereas possibility 
(ii) is realized for <f> > 4>d- In both cases, we are left after this stage with one monopole 
with left tail and E L = E^(I), and E R = E&\l) [case (i)], or E R = E^(I) [case 
(ii)], moving toward x = L. The current density is determined from the following 
equation, analogous to ( |4.13| ): 

dJ _ (£( fc > - £«) 2 v[ U k 
ds 



(4.46) 



where k = 2 in case (i) or k — 
with the initial condition I(sd 
in case (ii). 

In case (i), I'{sd) > if 



3 in case (ii). This equation has to be solved for s > Sd 
) in case (i), or for s > s m with initial data I(s m ) = v m 



(4.47) 



1< 6 < 



£7W v[+EW \v' 2 \ 



For the velocity curve we use, the bias interval (4.47) is very narrow, and the corre- 
sponding values of I(sd) are very close to 1. We may have a small-amplitude current 
oscillation corresponding to a supercritical Hopf bifurcation. When the bias is larger, 
I'(sd) < 0; and I decreases until / = v m is reached and Er becomes E^ 3 '(I). From 
there onwards, / increases and the situation is the same as in case (ii) after the fast 
intermediate stage. We then have a current oscillation of amplitude approximately 
given by 1 — v m . The transition from small to large- amplitude oscillation is extremely 
sharp, which accounts for the difficulty in observing it in simulations or in real lab- 
oratory experiments. In both cases, after some time the current density increases 
until the value Ip, the monopole flips its tail to the right, and we have completed 
the asymptotic description of one period of the self-sustained oscillation. Comparison 
between our asymptotic solution and a numerical simulation is shown in Fig. 5. 
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5. Birth of the new shock wave. In this Section we need to solve the char- 
acteristic equations given an approximate expression for the current density, e.g. 
Eq. ( 4. 28j ) , and then calculate the shock formation time and its earliest position. 
These quantities can be used as initial values for the following stage of the oscillation 
as described in the previous Section. 

There are two possibilities: 
• J(t s ) ~ 1 at the shock formation time. Then the current density may be 



approximated by Eq . fl4.28| ) during the process of shock formation. 
J(t s ) < 1, and Eq. (4.28) ceases to hold during the last part of the process 
of shock formation. After the current decreases below I = 1, we again have 
a quasisteady stage with 



(5.1) E L ~E^(I), E R ~EW(I), Y 



E^){I)- EM (I)' 



and 



(5.2) 



dJ 
ds 



(£(3) _ £(2))2 



■ Wo 



EW -(/>+(<]) -EW)v' 2 /v' 3 



<0, 



Clearly I(s) and Er ~ decrease while E- ~ increases as the time 
elapses. Eq. (5.2) approximates the current density from the time it has 
decreased below 1 until the shock formation time. If no shock wave is formed 
as approaches the value cf> (which occurs for c sufficiently close to 0), 
Y ~ 1 according to ( |5.1| ), and the old shock wave exits at the receiving end of 
the SL. At this time the stationary uniform field profile E = <f> is reached and 
maintained. These observations are in excellent agreement with the results 
of the numerical simulations. 



For a wide interval of c, Eq. ( 4.28 ) describes the current density during the full 
shock formation stage. We will see that the shock formation time decreases as c 
increases for small c, it reaches a minimum and then increases rapidly with c. Actually, 
Eq. ( |3.16| ) implies that the integral of v'(E(t;r s )) has to be negative, and therefore 
that the field on the characteristic with r = t s needs to take values on the second 
branch of the v(E) curve. Typically the field Eq(t) increases with r, it reaches a 
maximum, E™ ax , at T max and then decreases to E^(I(t)). We have to distinguish 
two cases depending on whether E™ ax surpasses 1, which happens for enough small c. 
In fact, the critical value, c, for which this maximum equals 1 occurs when the current 
density is maximum [observe that Eq. (4.28) does not depend on c\. As I max ~ 1 + ea, 
then it follows from ( [3.9[ ) that c = ea. 

Let us consider the case c < c in which this maximum value surpasses 1. Then 



the smallest time t s for which ( [3. 16 ) is satisfied, is the one corresponding to the 
characteristic issuing from Eq(t s ) — 1. For smaller t s there is initially a positive 
contribution to the integral, leading to higher values of t s . For larger r s the field along 
the characteristic quickly falls to E ^> (J) leading again to a positive contribution to 
the integral and therefore equation ( 3.1 6| ) is not satisfied. For E m > E(t;r) > 1 the 
integral of the characteristic equation along the direction of increasing t is unstable. 
Thus, we carry out the integration of the characteristic in the direction of decreasing 
t, starting at E^ (I (t s )) and ending at t such that E(t;r s ) = 1. The value of t s 
for which equation ( 3.16 ) is satisfied gives the shock formation time and therefore 
the shock position. In this case, t he sh ock formation time decreases with increasing 
c, because the right hand side of (|3.16|) increases. We have solved this problem for 
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L = 50, c = 1CT 4 obtaining t s - t M = 7.08, x s = 7.66 and I(t s ) = 0.93. This value 
agrees well with the result of direct numerical simulations of the model. 

If c > c, the smallest time t s for which ( 3.16] ) is satisfied, is the one corresponding 
to the characteristic issuing from Eq icix (t s — T max ). The characteristic corresponding 
to this value, r s , starts by giving a positive contribution to the integral until E(t;r s ) 
reaches 1. Therefore as c increases this positive contribution increases, and the shock 
formation time increases, which has been verified by direct numerical simulation. 

Let us now find approximations to Eqs. ( 3.16] ) - (3. IS) with the help of our previ- 
ous work on intermediate scales instead of integrating numerically the characteristic 
equations. Clearly we find s s — sm — 0(e 1 ^ 2 ), and y s — 0(e 1 ^ 2 ). While the old shock 



wave corresponds to a monopole with a right tail moving with velocity (4.7), the new 
shock wave corresponds to a monopole with a left tail, moving with velocity (4.6) (see 
Figure 3b). 

In order to estimate the shock formation time and position, we need to calculate 
the solution of ( ft,16| ) and (3.17) approximately. This might be quite involved, as 
the current density (needed to calculate the electric fields appearing on the formulas) 
depends on the field profiles, which should be calculated by the method of character- 
istics. Fortunately in the asymptotic limit e — > 0, the new shock wave is born during 
a time interval |s — sm\ = 0(y/e), at a distance y s = 0(*/e) from y = 0. The current 



density is approximately given by the parabola ( [4.281) , and both the equations for the 
field at the injecting contact and the characteristic equations are simpler. We shall 
give the approx imat e form of th ese equations and then insert their solutions in the 
exact formulas ( 3.16j ) and ( |3.17 ), whose solution yields the shock formation time. 
Let us start finding the approximate equation for the electric field at y — 0, (3.£). 



Durin g the quasistationary first stage of the oscillation described by Eqs. ( 4.11 ) to 
( 4.13 ), the field at the SL boundary is given by the quasi-steady form of ( |3.9| ), 



(5.3) 



E (s) ~ EW 



1(a) 



until s = s\ . During the faster stage about sm ch aracterized by the scaling ( 4.1E 
( 4.20 ), (3.9) may be approximated by using (fL27|) for I and E - 1 = e 1 / 2 E : 



(5.4) 



dE 



a — cL + 



K(i)|(i + C ) 



El 



a 2 s 2 



ds " 2 V u l + c, 

We can rewrite this equation so as to eliminate all parameters but one: 

de 



(5.5) 
(5.6) 

(5.7) 
(5.8) 



da 



[i + el - a 2 , 



fi = (l - — ) VT+c, 

ae 

e o — 7= -fro , 



2a 



(KqjH^q + c) 1 / 4 . 

V2 



This equation should be solved with the initial condition 

' 1 — /te N 



(5.9) 



E (s Kl ) - 



l + c 
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which follows from (5.3) when the current I(s K i) = 1 — ne (for s in the overlapping 
region) is u sed [see ( 4.30 )]. If c = O(e), we may solve (5.9) to leading order within 
the scaling (O)-(p): 



(5.10) 



('() 



/ n T\ K + cle 

(-y/l + K/a) = 



A phase plane study of Eq. (5.5) indicates that E increases, reaches a maximum 
and then decreases. There are several cases worth consideration according to the 
values of the paramet er /i (Figure 6). If c < ae, Eq surpasses 1 (Fig. 6a). In this 
case we ne ed t o solve (5.5) in order to find the initial conditi on f or the characteristic 
equations (5.5) with \i = 1. For smaller /i's, we may still use (|5.3| ) to approximate Eq. 
In either case the shock formation time will be obtained by solving the characteristic 
equations with an initial condition given by the solution of Eq. (5.5) or by (5.3). 

We now find the approximate form of the characteristic equations. In the scaling 
( |4.18D - J4.2CD they become: 



(5.11) 



dE 



|«"(1)| 



(E 2 



2 »2\ 

a s 



ds 2 

By redefining the variables in this equation, we can rewrite it as 

de\ 
da 



(5.12) 
(5.13) 
(5.14) 



1 



ei 



""(1)1 
2a 



En 



a\v"(l)\ 



Notice that a = (1 +c)- 1 / 4 <7 by (5.£). Whc n (7 takes on this value, the field on the 
characteristic is equal to the boundary field Eq. In scaled variables, the characteristic 
is parametrized by the value a such that when a = (1 + c) -1 / 4 <j, y = and 



(5.15) 



ei 



(l + c) 1 ^' (l + c )l/4 

The general solution of ( |5.12 ) is 



= (1 



n-3/4 



eo(o-). 



(5.16) 



ei(<r; a T ) 



L„ exp(-cr 2 + s 2 ) ds 



for a > Co (ei — > — oo as a 
coming from the initial condition (5.15). 

Inserting ( pU6| ) in (§J6|) and pl7D , we find 



where cto is a function of er r = a j(\ + c) 1 / 4 



(5.17) 



In 



e Ms 2 ~ ° 2 S ) ds 

1 \ 



Cy/T 



(5.18) 



L T exp(-cr2 + s 2) ; s y 

n 2/(l+c) 



dr 



Jan 



ds . 
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Here a s and a T correspond to t s and r s in Eq. (3.16) - (3.17), respectively. Notice 
that (To is a function of a T = <r/(l + c) 1 / 4 given by ( 5.15 ). We have solved numerically 
these equations for a fixed e and different values of c, obtaining a result consistent 
with that explained earlier: the shock formation time decreases as c increases for small 
c, it reaches a minimum and then increases rapidly with c. 

6. Death of the monopole. We only need to consider case (ii) separately, for 
case (i) does not require the introduction of different scales and was already analyzed 
in Section ^. In case (ii), an intermediate stage with the same time scales as those for 
the shock birth describes the death of the old shock wave. We make the ansatz 



(6.1) 
(6.2) 
(6.3) 
(6.4) 



I ~ v m + el, 
Em,r — E m ~ e 1 ^ 2 £ ; M,fl, 
E L -Eg>~eE L , E^=E^(v m ), 



a/2 



wher e s m i s the t ime at which the current reaches its minimum value. Equations 
B , ( 14361 ), and ([1371) become 



(6.5) 

(6.6) 

(6.7) 

(6.8) 
(6.9) 



dE M 
ds 

dE R 
ds 



— - f 1 / 2 II 
as 

dY 
ds 

= i-''- 
= i - 



+ G{e) 

+ 0(e) 

0(e^ 2 ) 

El + Oie 1 ' 2 ) 
i-v'{E$)E L =0{e^) 



e 1 ' 2 v 



E M 



where U m = U(E m ). Eqs. (|4.38| ) and ( [4.39|) still hold, but the approximation ( |4.40| ) 
does not. The solution of (|6. 



(6.10) 



^n(s) - ^n(Sm) + e 1 ' 2 U m s + 0(e). 



We now insert ( |4.33| ), Q6.2|) , ( |6.3| >, and ( p.lCD into (Mq ) and ( [4.391) . The result is 

E m — 4> 



(1) 



(6.11) 

where we have defined 
(6.12) 



1/2 [E R 5 -{8-1) E M ] (</> - E,, 



(E„ 



E 



(1)^2 



E,, 



Him 



E 



(i) 



(1-Y d 



Here Yd i s the position of the old shock wave when the current reaches v m , given by 
Eq. ([1.43]) with / = v m . Notice that 5 > 0; 6 — implies that the monopole Y arrives 
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at x = L exactly when the current density becomes v m . From the bias condition 
( |6.11 ), we find the position of the new shock at s = s m , and the relation 



(6.13) 
(6.14) 



bs 



E R S-(5 



DE 



M ■ 

(1)\2 



E { ') 



E 



(i) 



Notice that Em < < Er. Then ( 6.1 3| ) implies that 8 < 1, for otherwise we could 
not have s < 0. We now need to solve Equations (6.7), ( |3.S| ) and (6.13) for the 
unknowns Em, Er and I. Notice that ( |6 .13 ) also holds at the end of the previous 
quasistationary stage where two monopoles coexist. There Eqs. ( |6.7| ) and ([3^) with 
zero left side also hold, which gives us the following expression for the current at the 
end of the previous quasistationary stage: 



(6.15) 



V 

j _m 

2 



bs 



28- 1 



for < 8 < 1/2. 

We now find the unknowns Em, Er and /. From (6.7) and ( 5.13] ) we can find I 
and Em as functions of Er, which inserted into (3.E) yield 



(6.16) 



cLEr 



= b 



{Er — bs) [{l-26)E R + bs}, 



ds 2(1 - 5) 

We can rewrite this equation so as to eliminate all parameters but S 

df 



(6.17) 
(6.18) 



drj 



1 



/ = 



2(1 - 26) 



1/2 



{f-rj) [(1-26) f + V ], 

1/2 



E M , 



V = 



2(1-26) 



A phase plane study of (6.17) indicates that all trajectories tend to / = rj as r\ — > 
oo, which matches with the quasistationary stage that occurs after the disappearance 
of the old shock wave at X d (see below). This is illustrated by Figure 6d and by the 
exact expression, 



(6.19) 



f(v) = V 



1 



(1-25) flexp[(-^ + s^)6]ds 



where rjo is a constant. Of all these trajectories, the ones that match the adiabatic 
stage before the old shock dies are those that follow / = —77/(1 — 26) and experience a 
fast transition to / = rj. When / = —77/(1 — 25), the current matches ( |6.15 ). Assume 
that the fast transition takes place at 77 = r)\ 3> 1. An approximate expression for 
f(rf) may be found letti ng / — —frji/(l — 28), and 77 = 771 + 77', with 1 -C 77' <C r/i. 
Inserting this ansatz in ( |6.19 ), we find 



(6.20) 



df_ 
drf 



-r?! (f + l-26)(l-f)+0 



whose solution yields 
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where z is a constant. After the fast transition (7/ — > oo), both Em and Er become 
6s, the old shock dies and the current becomes 

(6.22) I ^ Vm + eb+'^b 2 (s-s m ) 2 . 

It is interesting to find out when the current / reaches the value v m . In terms of 
the variables rj and /, / may be written as 



(6.23) I = 



1 



Notice that I = if / = —rj/6. Then the current reaches / = v m (I = 0) if the slope 
of the line / = — rj/6 is larger than the slope of f{rj) as ry — > — oo, i.e. 1/8 > 1/(1 — 25). 
This implies < 6 < 1/3, which occurs for biases <fid < <j> < <j>s- When the bias is 
larger, 1/3 < 6 < 1/2, and the minimum of the current is larger than v m . 

7. Concluding remarks. We have developed an asymptotic theory for the con- 
tinuum limit of a discrete drift model of doped SLs. Our theory is in excellent quan- 
titative agreement with numerical simulations of the discrete model for long enough 
SLs. Our asymptotic analysis of the continuum equations indicates that they have 
a stable time-periodic solution for an appropriate bias. That this is indeed so, re- 
mains an important open mathematical problem. The time-periodic oscillation is due 
to repeated creation and propagation of traveling monopole wavefronts, formed by a 
shock wave and an attached tail region moving at the same velocity as the shock. The 
study of the shock-wave dynamics allows us to reduce the complexity of the problem 
by decomposing it into "coherent" structures and analyzing their simpler reduced dy- 
namics. This decomposition should in principle yield a full explanation of the classical 
chaos observed numerically when an appropriate harmonic forcing is added to the bias 

The present phenomena are related to the well-known Gunn instability in bulk 
semiconductors ]26|, ^4], [l3|, |[ fjj . The main mathematical difference is that the travel- 
ing waves responsible for the Gunn effect arc solitary waves of the electric field (charge 
dipolcs), while in the present case they are monotone wavefronts of the electric field 
(charge monopoles). We remark that the first numerical observation of a "Gunn- 
like" instability (in a drift-diffusion model with different boundary conditions from 
ours) mediated by monopoles is due to H. Kroemer, |l9|. An incomplete asymptotic 
study of the Kroemer drift-diffusion model valid for particular electron velocities was 
performed in [|l3). The fact that many different models (e.g., "slow" Gunn effect 
in ultrapure Germanium H, 0, 0) present the Gunn instability poses the problem 
of characterizing the model features that produce the Gunn effect, as precisely as 
possible. 

Other open problems presently being considered include explaining how well the 
asymptotic solution approximates the solution of the discrete equations and how well 
these results compare with experiments jl6[ |2Cj. To this end, it is important to 
extend the asymptotic analysis to a more complete model of SLs under photoexcitation 
|6|. Finally, a careful derivation of the discrete drift model from a fully quantum- 
mechanical setting would be most desirable. None of the presently known derivations 
is satisfactory. 

Acknowledgments. We thank J. Galan, H. T. Grahn, J. Kastrup, S.-H. Kwok, 
R. Merlin and A. Wacker for fruitful discussions and collaboration on related topics. 
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Appendix A. Brief Derivation of the Continuum Model from the Dis- 
crete Model. 

A SL is a periodic succession of alternating long (for simplicity) slabs of two 
semiconductors (GaAs and AlAs in jl6[ |2(J). Since the energy bandgap is different 
for the two semiconductors, the SL may be thought of as a periodic succession of 
potential "barriers" and "valleys". Along the growth direction of a finite SL, we find 
N periods formed by a barrier and a valley with a typical length of tens of nanometers. 
The lateral dimension may be one or ten thousand times larger, so that transport 
phenomena along the growth dimension may be considered to be one-dimensional. In 
the simplest case of a n-doped N-period SL without laser illumination, the discrete 
drift model consists of the following system of equations M, ff : 



(A.l) E j -E j _ 1 = ^{h j -N D ), 

dE 

(A.2) e —J- + e v(Ej) hj = J, 



dt 



N ^ J Nl 

3 = 1 



where j = 1, . . . , N. In this model Eq. (Al) and Eq. ( |A.2| ) are, respectively, the 



one-dimensional Poisson equation (averaged over one SL period) and Ampere's law 
for the electric field Ej(i) and electron density hj(t) at the site j, and for the total 
current density J(t). In these equations, the positive constants e, No, I and e are the 
average permittivity, average donor concentration, SL period and the absolute value 
of the charge of the electron. v(E) is an effective electron velocity to be specified 



later. Equation (A. 3) establishes that the average electric field is given by the dc 
voltage bias <!>. Notice that there are 2N+2 unknowns: Eq, E±, . . . , En, hi, ... , un, J 
and 2N+1 equations so that we need to specify one boundary condition for Eq plus an 
appropriate initial profile Ej(Q). The boundary condition for E (the average electric 
field before the SL) can be fixed by specifying the electron density at the first site, hi, 



according to (A.l). In typical experiments the region before the SL has an excess of 
electrons due to a stronger n-doping there than in the SL [|l6| ^(J . Thus it is plausible 
assuming than there is an excess number of electrons at the first SL period measured 
by a dimcnsionless parameter c: 

(A.4) El (})-E {t) = ---^. 

e 

c has to be quite small because it is known that a steady uniform-electric-ficld profile 
is observed at low laser illumination in undoped SL |2(], |j. If we adopt ( |A.4| ) with 
< c <C 1 as our boundary condition, a steady almost uniform electric field profile is 



clearly a possible solution of Eqs. (A.l) to (^ 

The only function not specified so far is the electron velocity v(E). This function 
is used to model sequential resonant tunneling (SRT) |(| in the SL, which is the main 
charge transport mechanism for the high electric fields. The velocity curve can be 
derived from experiments J^] or from simple one-dimensional quantum-mechanical 
calculations of resonant tunneling, as was done by Prengel et al |25| . For our analysis, 
the only important characteristic of the velocity curve is that it has to have a local 
maximum and a local minimum and therefore a region in which the velocity decreases 
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with increasing field. In our analysis we will use the curve depicted in the inset of 
Fig. la. 



It is convenient to render the equations (A.1)-(A.4) dimensionless by adopting as 
the units of electric field and velocity the values at the first positive maximum of the 
velocity curve, v(E), Em and vm (about 10 5 V/cm and 427 cm/s, respectively for the 
sample of Ref. Gq]). We set 



Ej 



Em 



(A.5) 



t = 



V_MVt_ 

I 



I = 



J 



eN D v M 
I 

N Em V 



where the dimensionless parameter v, defined as 

N D el 



(A.6) 



e Em 



is about 0.1 for the SL used in the experiments Jig ], Then the dimensionless equations 
of the model are: 



(A.7) 
(A.8) 

(A.9) 
(A.10) 



dE j 

-±+v{E J )n ] =E 



N 



Ei - E Q = v c; 



N 



c < 1. 



Here v(E) — v(E)/vm, and 

whereas v, c and N are fixed for each SL. Equations ( A.7 ) and (A 



is a dimensionless control parameter (the dc bias), 

are to be solved 



with initial conditions for the fields, Ej(0), compatible with the bias ( |A.9| ) and the 
boundary condition ( A.lCj ). The initial conditions for the electron density, ^(0), then 
follow from (A.7). 

In the continuum limit v —> 0, N — > oo so that 



(A.ll) 



= jv€ [0,L], L = Nv»l, 



the system of equations (A.7)-(A.10) becomes (|l.l|)-(|L~3|) after the electron density 
n(x,t) is eliminated by means of Poisson's equation. The hyperbolic equation (1.1) 
develops shock waves from initial data in finite time J2|, |24| . To determine their velocity 
and entropy condition we have to go back to the discrete model ( A.7)-( A.1C| ) B. The 
shock wave is a field profile that moves rigidly with an average speed V(E+, E-), so 
that Ej(t) = E(J-j s (t)) with E(-oo) = E-, E(+oo) = E + and j s (t + At) - j s (t) ~ 
V(E+, E-) At/v, for large enough At. j s (t) is the QW where the profile rij(t) reaches 
its maximum value at time t and X(t) — v j s (t) corresponds to the shock position in 
the continuum limit. Then Ej(t + At) — Ej(t) is approximately given by the distance 
js (t + At) — j s (t) that the shock has advanced during At times the field difference 
— (Ej — Ej-i) at some intermediate time in (t, t + At). Thus we have 



(A.12) 
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If we now sum (Ej — Ej_x) over the shock region, we obtain E^ 



hand rij — 0[y v ) ^> 1 and / = 0(1) in the shock region, so that ( A. 7 ) 
( |A.12D yield 



E_ . On the other 
) and 



E, - E- 



j=—oo 



( E 3 - E 3-l) 



(A.13) 



V(E+,E_ 



E 

j=-oo 



dEi 



v(Ej) dt 



E 

j=—oo 



Ea - E 



3-1 



In the continuum limit this becomes the equal area rule (1.4) 



E _ \v(E) V(E+,E_ 



dE = 0. 



To keep rij > inside the shock we must add the restrictions 

v{EJ) > V(E+,E-) > v{E + ) 

(and hence V = v(E) at only one point x inside the shock). Likewise, no shock with 
E+ < E- arises fro m re alistic initial conditions with rij > Jl3| . Typ ical solutions of 
the equal area rule (1.4) compatible with the entropy condition (1.5) are depicted in 
Fig. 2. Notice that we can have admissible values of E- and E+ belonging to equal 
or different branches of l/v(E). However for £L and E + to be on the same branch 
of l/v(E), they must belong to the second branch, with v'(E) < 0. 
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FIGURE CAPTIONS 

Figure 1. (a) Time evolution of the electric field profile on the SL using the velocity 
curve shown in the inset, (b) Charge density profiles, dE(x, t)/dx, showing the loca- 
tion of the IL for different times. The total current density versus time is shown in 
the leftmost inset, in which we have marked the times corresponding to the profiles 
depicted in part (a). The rightmost inset shows clearly a monopole with a right tail. 
Figure 2. Equal-area rule relating the electric field values behind, and in front, 
E + , of the shock wave with its velocity V(E + , E-). 

Figure 3. The time-periodic solution in a nutshell, (a) Motion of one shock during 
one period on a cylindrical surface whose axis is space. The diameter of the cylinder 
basis is 1 < 1 // < l/v m . (b) to (f) highlight of the equal- area conditions for the shock 
waves appearing in the electric field profile. They correspond to the different stages 
marked in Figure lb. 

Figure 4- Domain of the function E + = !F(E^, 1/V). The projection of E + (t) on this 
diagram for each time represents the evolution of the shock wave during its entire life. 
Figure 5. Comparison between our asymptotic solution and a numerical simulation 
of the full model for <j> = I .25, e = 0.02, and c = O.Of . 

Figure 6. Phase plane study of Eq. (5J3) corresponding to the field at x = during 
the stage of shock birth, (a) /i = f , which corresponds to the electric field on the 
characteristics that emanate from x = 0. (b) /i = 0.5. (c) [i — —0.5. (d) P hase plane 
describing the death of a shock wave when / v rn . It corresponds to Eq. (6.17) with 
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